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f^ ' Quantum Monte Carlo methods have proven to predict atomic and bulk properties of light and 

flj , non-light elements with high accuracy. Here we report on the first variational quantum Monte Carlo 

^O ' (VMC) calculations for solid surfaces. Taking the boundary condition for the simulation from a finite 

layer geometry, the Hamiltonian, including a nonlocal pseudopotential, is cast in a layer resolved 

form and evaluated with a two-dimensional Ewald summation technique. The exact cancellation 

of all Jellium contributions to the Hamiltonian is ensured. The many-body trial wave function 

consists of a Slater determinant with parameterized localized orbitals and a Jastrow factor with a 

common two-body term plus a new confinement term representing further variational freedom to 

C/3 ' take into account the existence of the surface. We present results for the ideal (110) surface of 

Galliumarsenide for different system sizes. With the optimized trial wave function, we determine 

some properties related to a solid surface to illustrate that VMC techniques provide standard results 

under full inclusion of many-body effects at solid surfaces. 
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I. INTRODUCTION 



o, 

O , The theoretical description of solid surfaces bears fundamental aspects as well as technologically important ap- 

plications. Most of the theoretical approaches are based on common local density functional theory (DFT-LDA) or 
its improvements (GGA). Despite its success a systematic development is hampered by the inherent approximations 
^ . made, especially for highly correlated or inhomogeneous systems. Therefore, alternative approaches are desirable to 
C — ' test and further develop the quahty of the exchange-correlation functional. It is a long tradition to use Quantum Monte 
'"\l" , Carlo (QMC) techniques for the homogeneous electron gasld not only to describe the system itself, but also to appj^i 
the results as an approximation to inhomogeneous systems. Especially through the use of nonlocal pseudopotentiala3'0 
within variational quantum Monte Carlo (VMC), the extension of QMC methods from homogeneous systems to soL 
^^ • of non-light elements is state flttJie art with a variety of aspects such as cohesive properties of elementala or bipolaii 
f~^ I solids, or excitation energies. LfcliJ 
■^ . Here we propose a scheme to apply the VMC method to solid surfaces. Being safely founded on Ritz's variational 

Cd ' principle, all approximations made in the trial wave function are controlled. An unknown functional as in DFT or a 
H , perturbative approach like a diagrammatic many-body expansion as in GW is thereby avoided. 
I ■ The paper is organized as follows. In Sec. ||, we describe the method with a detailed discussion of the boundary 

' ^ ' conditions, the Hamiltonian and the trial wave function. In Sec. [II, we present results for the optimization of the 
« 
O 
O 



confinement term for two different system sizes keeping other surface specific parameters in the one-body orbitals 
fixed. The dependence of the total energy on the coordinates of a surface ion is examined to test the sensitivity of 
the method to surface specific information. 



lJ ■ II. METHOD 

P.. 

A. Stochastic Integration 

In VMC the expectation value of the Hamiltonian of a quantum system in the state ^> 



J\^x{Y)\ dY *A(r) 



is calculated by stochastic integration over the coordinates Y of the system, where A denotes a set aLparameters 
representing the variational freedom in the trial wave function. In practice the Metropolis algorithir£3 is used to 
guarantee importance sampling during the random walk so that the energy can be estimated safely as the mean of 
the local energy with purely statistical error 



Ritz's variational principle ensures that the minimum Ex-^ of Ex in the parameter space is an upper bound for the 
true ground state energy Eq. 

B. Boundary conditions for a solid surface 

The central problem in a QMC simulation of a solid surface is the choice of the boundary conditions along the 
surface normal which is chosen here to be the z-axis of the coordinate system. Often a supercell geometry is used to 
model the semi-infinite crystal: the simulation cell contains a slab of layers surrounded by vacuum on both sides in 
the z-direction, and periodic boundary conditions are applied to all three directions. The use of a localized basis eases 
the choice of suitable boundary conditions. We assume periodic boundary conditions only in the surface plane and 
take a symmetric slab of layers with the normalizability of the wave function as the only boundary condition in the 
z-directi on. T his layer resolved formulation c an be generalized to the semi-infinite crystal as well. The Hamiltonian 



(see Sec. |IIC| ) and the wave function (see Sec. 11 D) are both affected by the choice of the boundary conditions. First, 
we describe the adaption of the geometry of the simulation cell to the (110) surface of GaAs. 

With the z-axis along the [110]-direction and the x,y-axes in the (llO)-plane, the following set of vectors spans 
the lattice of GaAs: ai = a (4=, 0,0), 02 = a (0,1,0), 03 = '^{^t^t^)^ the basis consisting of tas = (0,0,0) and 

TGa — 0.(^1 jiO)i where a is the lattice constant of cubic GaAs. The two-dimensional lattice {Ru} of the simulation 

cell is chosen as the span of &i = NX • ai, 62 = NY • a2. Denoting the number of lattice planes inside the simulation 
cell by NZ, we completely characterize the cell by the triple (NX, NY, NZ) containing N = 8 • NX • NY • NZ electrons 
and K = 2 ■ NX • NY • NZ ions. NZ must be large enough to prevent artificial interactions between the two surfaces. 
Formally the size of the simulation cell in z is infinite because no further assumptions about the z-direction are made. 
That is why no interaction between the two surfaces through the vacuum can occur as in the usual supercell-slab 
geometry. 

To model the infinite solid in this geometry, taking as the third primitive vector 63 — NZ • (a3)_L — NZ • a • (0, 0, -y=) 
gives an orthorhombic simulation cell Q with the constraint NZ to be even. 

C. Hamiltonian 

The Hamiltonian of the considered system, here electrons of a solid or of a solid surface in the Born-Oppcnheimcr 
approximation, consists of the kinetic energy operator T of the electrons and the potential energy operator V of all 
the Coulomb interactions in the system. For solids of non-light elements the inner shell electrons and the nucleus 
can be replaced by a pseudo-ion with effective charge Z and a nonlocal electron-ion pseudopotential for the valence 
electrons.l3 The potential energy operator reads 

V^Vc + Vhp + V^'p , (3) 

where Vc is the operator of the remaining Coulomb interactions and T/^™-and Vpp are the local resp. nonlocal part of 
the pseudopotential. The nonlocal part is evaluated in a semi-local form.Ll The resulting Coulomb energy like operator 
collects together the constant ion-ion-energy, some part of the electron-ion interaction and the full electron-electron 
interaction 

Vc = Vr + Ve^ + Vee , (4) 

thereby including the contributions of the Z/r-asymptotic of Vpp among Vd- The remaining pseudopotential part is 
short ranged, and only Vc is affected by the boundary conditions. To account for the surface plane parallel periodic 



boundary conditions, a super-lattice with the vectors {Ru} is introduced, which adds to an arbitrary position of 
charge its periodic counterparts. For clarity and later reference, Vc explicitly reads 



2Vc := E 
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where a runs over the index set of gallium and arsenic ions and i over the electrons in the simulation cell. As in the 
three-dimensional case each of these two-dimensional lattice summations over u diverges for every r. 

One possibility to overcome this is to introduce a neutralizing background charge in every interaction. To ensure 
that all these jellium contributions do not artificially change the total Coulomb energy, the background charges are 
forced to cancel out locally, at every point, by adding the following terms to Vc^ which vanish due to charge neutrality 
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where R = (i?|| ,Rz)- Splitting up the r. h. s. of eq. (^ into two terms, these are fed into the first resp. second sum over 
i of eq. (|^). In the same way the parts of eq. (0) are distributed over the sums over a. Now all interactions of eq. (J^) 
converge and can be calculated separately. We define the following two-dimensional surface lattice potential v" 



v%r) := y( ^ + r{r)) - i- 



(8) 



with 
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where J7q is the area of the surface unit cell, and zq chosen so that J7 = J7g x [—^,+^]. By generalizing ideas of 
Ref. |lj we evaluate w^ by Ewald summation technique in the form 
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The result for Vc for the simulation of the surface system reads 



V^c4E(-^(-.)) + ^E^ (11) 

^E(^"-^(-»))-E|i 

I. a I, a 

+ ^E(^"^"-^(--')) + ^eS^ 

+ v:/({^.},{^a}), 

= V; + Kj , (12) 

where Vah — fa — fb a-nd z^, — {Ra)z- The Coulomb interaction is decomposed into a term V^ depending on the 
inter-particle difference vector only, which gives the main contribution to the total energy, and into a term Vj which 
collects up all remaining terms. As is shown in appendix A, these depend on z-coordinates only. Appendix B contains 
details of the numerical implementation. 

The everywhere bounded, relatively smooth surface lattice potential v'^ is used per tabulated interpolation with a 
Herman-Skilman mesh in z-direction. The tabulation is necessary to make practical calculations feasible. 

D. Wave function 

Our ansatz for the wave function ^ = $5 x vj/j consists of a product "^ s of Slater determinants and a Jastrow 
factor \E'j with a two-body term and a confinement term 

*(fio-i, . . .,fN(7N) = D^ri, . . . ,f«)D^(f«+i, . . . ,r7v) 

X exp (^ - ^ u{nj,aij) - Y^ Vcf{fi)j (13) 

i<j i 

with Gij denoting parallel or anti-parallel spins. The two-body term u is chosen in a periodic form 

u{f,a)^ A{^+v''{f)\ (l -€-"■' ^^"A . (14) 

A is a variational parameter and F{a) is, for given A, fixed by the cusp conditiorcj to compensate for the Coulomb 
singularity. 

The original formulation of Fahy, Wang, and Louieaa included a one-body function xi^i) in order to correct for the 
effect of the two-body term u of smoothing the electron charge density, and to restore the LDA charge density. This 
idea has been extended successfully to enlarge the variational freedom in the Jastrow factor keeping the one-body 
states in the Slater determinant unchanged.Q But one drawback of this approach is that the node surface of ^ is 
entirely determined by ^s so that e.g. further diffusion quantum Monte Carlo (DMC) calculations are restricted to 
the quality of the node surface of the underlying VMC wave function which in fact depends on the quality of the 
one-body orbitals. _ 

Our approach allows for a direct optimization of the one-body orbitals in the Slater determinantu viz. 

'l^o.,{f) = <Ptl{f) + P€l{f) , (15) 

<t>Z(f) = 7^^'^f ([r - i?^^]/e) + E ^;<'([^^- Rt'^yc.^i , 
&f) = 7°^0?^([r - i?f 1/cf ^) - E ^l^^'ii^- ^f I/O 

i 

for a zincblende system as GaAs. The index a here only runs over the set of arsenic ions in the simulation cell, /i 
over the four tetrahedral directions 17, originating from an arsenic ion, and i over the Cartesian coordinates (x, y, z). 
The index a denotes the index of the corresponding gallium ion for each bond and is fixed for given a and /i. These 
localized orbitals are parameterized hybrid bonds with seven variational parameters /3, 7, Q in the bulk case. Especially 
the contraction parameters C, allow for compensation of the Jastrow factor. The unmodified orbitals (/)*^^, 4>^^ are 4s- 
and 4p- pseudo wave functions coewsoonding to the nonlocal ab-initio LDA, generalized norm-conserving atomic 
pseudopotential following Hamann.tBEJ 



In the construction of the locahzed one-body orbitals for the arsenic ions in the outermost layer the corresponding 
gaUium bond partners are missing and new parameterized surface states have to be set. 

Our ansatz for the one-body surface states consists of a doubly occupied arsenic dangling bond orbital. It formally 



arises from the bulk orbitals of eq. ( 15) by setting (3 — if a denotes an arsenic ion in the outermost layer and /i — jidb 



denotes a tetrahedral direction towards the vacuum. Apart from 7?'^ , CP^ , and dp^ as three new surface specific 
composition and contraction parameters we allow for the variation of the direction of the dangling bond orbital by 
introducing two parameters, the angles "d and 0, determining the prefactor t* ^. To allow for effects of relaxation, i.e. 
the deviation of the ions from their ideal positions, in the construction of the hybrid orbitals the prefactors t' are 
chosen proportional to the vector connecting two neighboring ions. 

With this variational freedom one would expect that increasing the number of layers the total energy per atom 
comes asymptotically closer to the theoretical bulk value, because the surface to bulk fraction decreases. But, as will 



be stated in Sec. Ill B, the energy actually increases instead of decreasing with increasing system size. To compensate 
for the observed overshoot in smoothing the electron density over the whole range of the simulation cell by the Jastrow 
factor, we introduce the following confinement term in the Jastrow factor 

Vcf{T^ = Vcf{z) = - ■ Vcf [tanh(S'c/(|2:| - z^j - L^f)) + ij , (16) 

where Zcf is the z-distance of the outermost layer from the center. 

III. RESULTS AND DISCUSSION 

A. Bulk 

In order to test the numerical implementation of the layer resolved formulation of the Hamiltonian, of the geometry 
of the simulation cell, and of the construction of the trial wave function by comparison with independently obtained 
results, a bulk GaAs system is modeled by an orthorhombic cell called (2,2,8) containing 256 electrons. Usual three- 
dimensional boundary conditions are applied for the random walk, but two-dimensional boundary conditions for the 
Coulomb potential and a direct summation over the copies in z-direction, which can be truncated after a small number 
of copies because the potential is vanishing rapidly with increasing z-distance. The total energy is determined for 
different cubic lattice constants a keeping all other variational bulk parameters fixed. The resulting minimal energy 
per unit cell E^ = (—232.92 ± 0.05) eV/u.c. at an equilibrium lattice constant of Oq = (10.69 ± 0.05) a.u. is in very 
good agreement with our former results of E'^°^ = (—233.04 ± 0.08) eV/u.c. at an equilibrium lattice coastant of 



a, 



rcf 




(10.69 ± 0.1) a.u. obtained within the usual VMC scheme for solids for the same number of electrons.l3 



We conclude from this comparison that, with respect to numerical accuracy, the layer resolved implementation is 
suitable to tackle the many-body problem for solid surfaces. 

B. Solid Surface 

First, VMC calculations without confinement terms were performed. Only technically necessary changes due to the 
new form of the simulation cell, due to the prescribed two-dimensional periodicity concerning the Hamiltonian, the 
trial wave function and the random walk, and due to the lack of binding partners for ions at the outermost layers, 
were made. One would expect that with increasing number of layers the energy should converge to the asymptotic 
of the (theoretical) bulk energy E^. But for the surface system (2,2,8) an energy of E^ = (—154.13 ± 0.16) eV/u.c. 
is obtained. What is even worse, increasing the number of layers up to 12, the total energy increases monotonically 
with an average slope of 24 eV/layer. So a straightforward, brute force usage of the traditional VMC scheme seems 
not possible. 

Analyzing the different contributions to the total energy, one finds that with respect to the bulk case, a per unit 
cell decrease in the positive electron-electron interaction for larger systems is drastically over-compensated by an 
increase in the negative electron-ion interaction. This trend is even sharpened for bigger system sizes. Such drastic 
changes seem to be caused by electrons drifting extraordinarily far into the vacuum region. The actual existence of 
this mechanism is proved by observation of the z-resolved density. As one can see from the left part of Fig. (|l|) the 
electron density is smeared out into the vacuum regions. This behaviour comes from the known effect of the two-body 
part of the Jastrow factor of smoothing the electron density. In the bulk case, we have compensated for this effect 
through the introduction of contraction parameters in the localized orbitals of the Slater determinant. But here. 



having allowed the outermost, arsenic dangling bond orbitals to be contracted, we find no significant improvement on 
the relevant scale of 10 eV. 
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FIG. 1. z-resolved density with (right) and without (left) confinement term in the trial wave function 



Taking advantage of the flexibility of VMC as a variational procedure, we introduce an additional factor in the trial 
wave function allowing for a compensation of the Jastrow pressure on the scale of the whole surface system. Formally, 
this factor belongs to the Jastrow term in the form of a z-dependent confinement function Vcf- 

To examine the behaviour of the total energy with respect to the variation of the three confinement parameters, 
calculations for the system (2,2,8) with 256 electrons were performed yielding a minimal energy of i^lsg — (—223. 00 ± 
0.99) eV/u.c. for (K/, 5'^/, L^/) = (15.15 ± 1.64,0.137 ± 0.013,3.01 ± 0.77). This is an essential improvement in 
comparison to the former results. The optimization was done by a weighted regression where the total energy is 
written as a quadratic function of the variational parameters. 

In the region of parameter space where the global minimum is actually located, the information of one of the 
three parameters is redundant due to the presence of noise in form of the statistical uncertainty of the Monte Carlo 
data Ex. It is sufficient to fix one parameter at an approximately correct value. E.g. setting Lcf — 3.0 a. u. and 
scanning the remaining two-dimensional parameter space yields a fitted energy of -E256 — (—222.92 ±2. 16) eV/u.c. for 
{Vcf,Scf,Lcf) = (15.27 ±0.24, 0.139 ±0.006, 3.0). Finally, gathering statistics yields i^lgg = (-222.86 ±0.05) eV/u.c. 
We conclude that the actual uncertainty in energy is of the order of 100 meV/u.c. and that the trial wave function 
is essentially improved. Computed with the optimized trial wave function, now the z-resolved density as seen in the 
right part of Fig. iMj is much more localized to the solid proving that the above reasoning to explain the failure of 
the bulk derived trial wave function was correct. 

Whether the new confinement term is generally suitable also for larger systems, the trend of increasing energy with 
number of layers should be reversed. Therefore we performed calculations for the system (2,2,12) with 12 layers and 
384 electrons. Again a two-dimensional scan of the total energy landscape was performed yielding an fitted optimal 
energy of i^|84 = (-223.64 ± 2.09) eV/u.c. for (K/, S'c/, Lc/) = (32.68 ± 0.21,0.103 ± 0.001, 3.0). Repeated sampling 
with these values of parameters gives E^g^^ = (—223.52 ± 0.06) eV/u.c. 

So by increasing the system size a slight decrease in total energy is found which was originally to be expected from 
physical reasoning. Comparing the optimal values for the variational parameters underlines that the bigger systems 
needs a higher compensation of the Jastrow pressure which is reflected in the corresponding values of Vcf and Scf- 

Further lowering in total energy, now on the scale of 1 eV, has to be achieved by other variational degrees of freedom, 
e.g. by variation of parameters in the one-body orbitals in the Slater determinant. These results will be presented in 
a forthcoming publication. 

Here in the rest of this section we focus on the question, to which extent the proposed VMC procedure can 
resolve surface specific information. As a test case an in-layer displacement of the arsenic ions in the outermost 
layer is chosen. Is this statistically resolvable within the current stage of optimization of the trial wave function? To 
improve the fiexibility of the trial wave function with respect to ion displacement, the coefficients of the pi-orbitals 
are chosen so that the resulting hybrid bond orbitals are automatically concentrated in the bond direction keeping 
all other variational surface specific parameters fixed. In Fig. (^) the total energy per unit cell is plotted against the 
displacement of an arsenic ion perpendicular to the As-Ga chains of the ideal (110) surface. Despite the presence of 
statistical noise a minimum is clearly detectable. The equilibrium position is found very close to the ideal position 
slightly moved to the direction away from the As-Ga chains. So the determination of elements of the dynamical matrix 



via a parabolic fit is statistically relevant within this VMC procedure. 
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FIG. 2. Total energy in [eV/atom] of slab of 8 layers with outermost arsenic ions displaced within layer plotted against 
displacement in atomic units. Error bars denote the statistical uncertainty a{E\) of individual Monte-Carlo energies E\. Solid 
line shows the best fitting parabola. 

In conclusion, we proposed a scheme to describe solid surfaces in the frame of VMC techniques in order to use 
their fundamental aspects also in this field. Starting from the bulk derived, only slightly modified trial wave function, 
showed a failure with respect to finite size analysis. The reason was found in a large drift of the electron density into 
the vacuum. By introducing a new confinement term to the trial wave function much better results, with respect to the 
absolute value in total energy and with respect to the size dependence, were obtained. The sensitivity and relevance 
of the technique to surface structures could be proved showing that a surface analysis is feasible. Further systematic 
surface specific improvements in the trial wave function should allow to start from a well founded, variational basis 
to finally yield surface specific statements, as relaxation or reconstruction, and many-body influences therein. 
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APPENDIX A: EXPLICIT FORM OF Vj 

After insertion of eqs. (^) and (|^) into eq. (||) and regrouping, four terms like 
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occur, where 7 can either run through an electron index set or an ion index set, read j — a resp. 7 = z in cq. dq). 
The prime at the summation over 1/ means that the divergent term l/j/^fj'^'^l has to be skipped in the ^ = term of 
the sum if the outer indices 7 and 7' are equal and denote the same particle, however, the associated integral still 
being kept. 

Substituting R' := RY. + R and skipping the prime, the volume of integration changes from 17 = IIq x [—^ , +^'1 
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The next step is to substitute R = R' + f^' in the integrals and skipping the prime again. Now the integration runs 
over fi^o ■= (X]j/ ^t) X [~^ ~ ^1' ' +^ ~ •^7']- Furthermore, using the additivity of integration yields 
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Integration over fl'^ — f^oo means that contributions from ll^o count positively while f^oo contributes negatively. The 
first line of eq. (A4) depends only on the difference vector r^^i := r\ 



I ^ - r^i . So it is useful to define 
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Though rty and 7%^/ formally appear on the r. h. s. of eq. (A6), the actual dependence on r\\ drops out due to the 

infinite integration parallel to the surface. Changing the variable of integration to R' = R + r^i in eq. ( A6) leads to 
the result 
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with W^ := {J2i, ^t) X [~^ + -^7' ' +^ + -^7']- Physically this corresponds to a dipole like charge configuration with 
films of finite width. A co mpa ct solution of this problem of electrostatics is given in appendix B. Finally, inserting 
eqs. ( |A5| ) and (A?) in eq. (A_4), with additional prefactors Z^, Z^i to account for arbitrarily charged particles in the 
general case, 2T4fc can be expressed as 
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With these abbreviations the total Coulomb energy of eq. (H) can be written as 



^vi^ + wi^ + v^ + vi + v^ + vi + vd 

=■■ Vs + Vj. 



(A8) 
(A9) 

(AlO) 
(All) 
(A12) 



APPENDIX B: NUMERICAL EVALUATION OF Vj 



As defined in eq. (A9) 2Vj consists of terms like 



(Bl) 



Due to the symmetry of the problem (f-^ can be expressed by 
Before giving the explicit expression for f'-{ we define 
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Then the potential generated by a dipole like charge distribution of two oppositely charged films each of thickness 



dp with their inner borders separated by dg is given by 
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The potential is determined up to a constant Lp^. This is fixed by the convention that a charge at z = — oo feels zero 
potential leading to ip^ = Q. 
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